//不对齐
#include<iostream>
#include<arm_neon.h>
#include<sys/time.h>
#include<unistd.h>
using namespace std;
const int _NUM = 1024;
float matrix[_NUM][_NUM];

void Generate(float A[][_NUM], int size) {
	for (int i = 0; i < size; i++) {
		for (int j = 0; j < i; j++)
			A[i][j] = 0.0;
		A[i][i] = 1.0;
		for (int j = i + 1; j < size; j++)
			A[i][j] = rand() % 10;
	}
	for (int k = 0; k < size; k++)
		for (int i = k + 1; i < size; i++)
			for (int j = 0; j < size; j++)
				A[i][j] += A[k][j];
}

void Gauss_Neon(float A[][_NUM], int NUM) {
	for (int k = 0; k < NUM; k++) {
		float32x4_t v1 = vmovq_n_f32(A[k][k]);
		int j;
		for (j = k + 1; j + 4 <= NUM; j += 4) {
			//A[k][j] /= A[k][k];
			float32x4_t va = vld1q_f32(A[k] + j);
			va = vdivq_f32(va, v1);
			vst1q_f32(A[k] + j, va);
		}
		for (; j < NUM; j++)
			A[k][j] /= A[k][k];
		A[k][k] = 1.0;
		for (int i = k + 1; i < NUM; i++) {
			float32x4_t vaik = vmovq_n_f32(A[i][k]);
			for (j = k + 1; j + 4 <= NUM; j += 4) {
				//A[i][j] = A[i][j] - A[i][k] * A[k][j];
				float32x4_t vakj = vld1q_f32(A[k] + j);
				float32x4_t vaij = vld1q_f32(A[i] + j);
				float32x4_t vx = vmulq_f32(vaik, vakj);
				vaij = vsubq_f32(vaij, vx);
				vst1q_f32(A[i] + j, vaij);
			}
			for (; j < NUM; j++)
				A[i][j] = A[i][j] - A[i][k] * A[k][j];
			A[i][k] = 0;
		}
	}
}
int main() {
	Generate(matrix, _NUM);
	double sum_time = 0.0;
	struct timeval start;
	struct timeval end;
	int ecx = 10;
	while (ecx > 0) {
		Generate(matrix, _NUM);
		gettimeofday(&start, NULL);
		Gauss_Neon(matrix,_NUM);
		gettimeofday(&end, NULL);
		sum_time += 1000000 * (end.tv_sec - start.tv_sec) + end.tv_usec - start.tv_usec;
		ecx--;
	}
	cout << "The algorithm takes:" << sum_time / double(1000.0) / double(10.0) << "ms" << endl;
	return 0;
}
/*#include<iostream>
#include<arm_neon.h>
#include<sys/time.h>
#include<unistd.h>
using namespace std;
const int _NUM = 1024;
float* matrix;

void Generate(float A[][_NUM], int size) {
	for (int i = 0; i < size; i++) {
		for (int j = 0; j < i; j++)
			A[i][j] = 0.0;
		A[i][i] = 1.0;
		for (int j = i + 1; j < size; j++)
			A[i][j] = rand() % 10;
	}
	for (int k = 0; k < size; k++)
		for (int i = k + 1; i < size; i++)
			for (int j = 0; j < size; j++)
				A[i][j] += A[k][j];
}

void Gauss_Neon(float A[][_NUM], int NUM) {
	for (int k = 0; k < NUM; k++) {
		float32x4_t v1 = vmovq_n_f32(A[k][k]);
		int j;
		for (j = k + 1; j + 4 <= NUM; j += 4) {
			//A[k][j] /= A[k][k];
			float32x4_t va = vld1q_f32(A[k] + j);
			va = vdivq_f32(va, v1);
			vst1q_f32(A[k] + j, va);
		}
		for (; j < NUM; j++)
			A[k][j] /= A[k][k];
		A[k][k] = 1.0;
		for (int i = k + 1; i < NUM; i++) {
			float32x4_t vaik = vmovq_n_f32(A[i][k]);
			for (j = k + 1; j + 4 <= NUM; j += 4) {
				//A[i][j] = A[i][j] - A[i][k] * A[k][j];
				float32x4_t vakj = vld1q_f32(A[k] + j);
				float32x4_t vaij = vld1q_f32(A[i] + j);
				float32x4_t vx = vmulq_f32(vaik, vakj);
				vaij = vsubq_f32(vaij, vx);
				vst1q_f32(A[i] + j, vaij);
			}
			for (; j < NUM; j++)
				A[i][j] = A[i][j] - A[i][k] * A[k][j];
			A[i][k] = 0;
		}
	}
}
int main() {
	matrix = (float*)_aligned_malloc(_NUM * _NUM * sizeof(float), 32);
	//再将一维数组转为二维数组
	float(*matrix2D)[_NUM] = (float(*)[_NUM])matrix;
	Generate(matrix2D, _NUM);
	double sum_time = 0.0;
	struct timeval start;
	struct timeval end;
	int ecx = 10;
	while (ecx > 0) {
		Generate(matrix2D, _NUM);
		gettimeofday(&start, NULL);
		Gauss_Neon(matrix2D,_NUM);
		gettimeofday(&end, NULL);
		sum_time += 1000000 * (end.tv_sec - start.tv_sec) + end.tv_usec - start.tv_usec;
		ecx--;
	}
	cout << "The algorithm takes:" << sum_time / double(1000.0) / double(10.0) << "ms" << endl;
	_aligned_free(matrix);
	return 0;
}
*/